{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 7\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *\n",
    "\n",
    "from pandas import read_html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from the previous chapter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'data/World_population_estimates.html'\n",
    "tables = read_html(filename, header=0, index_col=0, decimal='M')\n",
    "table2 = tables[2]\n",
    "table2.columns = ['census', 'prb', 'un', 'maddison', \n",
    "                  'hyde', 'tanton', 'biraben', 'mj', \n",
    "                  'thomlinson', 'durand', 'clark']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "un = table2.un / 1e9\n",
    "un.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "census = table2.census / 1e9\n",
    "census.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_results(census, un, timeseries, title):\n",
    "    \"\"\"Plot the estimates and the model.\n",
    "    \n",
    "    census: TimeSeries of population estimates\n",
    "    un: TimeSeries of population estimates\n",
    "    timeseries: TimeSeries of simulation results\n",
    "    title: string\n",
    "    \"\"\"\n",
    "    plot(census, ':', label='US Census')\n",
    "    plot(un, '--', label='UN DESA')\n",
    "    plot(timeseries, color='gray', label='model')\n",
    "    \n",
    "    decorate(xlabel='Year', \n",
    "             ylabel='World population (billion)',\n",
    "             title=title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Simulate the system using any update function.\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that computes the population next year\n",
    "    \n",
    "    returns: TimeSeries\n",
    "    \"\"\"\n",
    "    results = TimeSeries()\n",
    "    results[system.t_0] = system.p_0\n",
    "    \n",
    "    for t in linrange(system.t_0, system.t_end):\n",
    "        results[t+1] = update_func(results[t], t, system)\n",
    "        \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quadratic growth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the implementation of the quadratic growth model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func_quad(pop, t, system):\n",
    "    \"\"\"Compute the population next year with a quadratic model.\n",
    "    \n",
    "    pop: current population\n",
    "    t: current year\n",
    "    system: system object containing parameters of the model\n",
    "    \n",
    "    returns: population next year\n",
    "    \"\"\"\n",
    "    net_growth = system.alpha * pop + system.beta * pop**2\n",
    "    return pop + net_growth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a `System` object with the parameters `alpha` and `beta`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_0 = get_first_label(census)\n",
    "t_end = get_last_label(census)\n",
    "p_0 = census[t_0]\n",
    "\n",
    "system = System(t_0=t_0, \n",
    "                t_end=t_end,\n",
    "                p_0=p_0,\n",
    "                alpha=0.025,\n",
    "                beta=-0.0018)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_simulation(system, update_func_quad)\n",
    "plot_results(census, un, results, 'Quadratic model')\n",
    "savefig('figs/chap07-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Can you find values for the parameters that make the model fit better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equilibrium\n",
    "\n",
    "To understand the quadratic model better, let's plot net growth as a function of population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop_array = linspace(0, 15, 100)\n",
    "net_growth_array = system.alpha * pop_array + system.beta * pop_array**2\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('whitegrid')\n",
    "\n",
    "plot(pop_array, net_growth_array)\n",
    "decorate(xlabel='Population (billions)',\n",
    "         ylabel='Net growth (billions)')\n",
    "\n",
    "sns.set_style('white')\n",
    "\n",
    "savefig('figs/chap07-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like.  Remember that the x axis is population now, not time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the growth rate passes through 0 when the population is a little less than 14 billion.\n",
    "\n",
    "In the book we found that the net growth is 0 when the population is $-\\alpha/\\beta$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "-system.alpha / system.beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the equilibrium the population tends toward."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sns` is a library called Seaborn which provides functions that control the appearance of plots.  In this case I want a grid to make it easier to estimate the population where the growth rate crosses through 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dysfunctions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When people first learn about functions, there are a few things they often find confusing.  In this section I present and explain some common problems with functions.\n",
    "\n",
    "As an example, suppose you want a function that takes a `System` object, with variables `alpha` and `beta`, as a parameter and computes the carrying capacity, `-alpha/beta`.  Here's a good solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity(system):\n",
    "    K = -system.alpha / system.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity(sys1)\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see all the ways that can go wrong.\n",
    "\n",
    "**Dysfunction #1:** Not using parameters.  In the following version, the function doesn't take any parameters; when `sys1` appears inside the function, it refers to the object we created outside the function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity():\n",
    "    K = -sys1.alpha / sys1.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity()\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This version actually works, but it is not as versatile as it could be.  If there are several `System` objects, this function can only work with one of them, and only if it is named `system`.\n",
    "\n",
    "**Dysfunction #2:** Clobbering the parameters.  When people first learn about parameters, they often write functions like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity(system):\n",
    "    system = System(alpha=0.025, beta=-0.0018)\n",
    "    K = -system.alpha / system.beta\n",
    "    return K\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity(sys1)\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we have a `System` object named `sys1` that gets passed as an argument to `carrying_capacity`.  But when the function runs, it ignores the argument and immediately replaces it with a new `System` object.  As a result, this function always returns the same value, no matter what argument is passed.\n",
    "\n",
    "When you write a function, you generally don't know what the values of the parameters will be.  Your job is to write a function that works for any valid values.  If you assign your own values to the parameters, you defeat the whole purpose of functions.\n",
    "\n",
    "\n",
    "**Dysfunction #3:** No return value.  Here's a version that computes the value of `K` but doesn't return it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity(system):\n",
    "    K = -system.alpha / system.beta\n",
    "    \n",
    "sys1 = System(alpha=0.025, beta=-0.0018)\n",
    "pop = carrying_capacity(sys1)\n",
    "print(pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function that doesn't have a return statement always returns a special value called `None`, so in this example the value of `pop` is `None`.  If you are debugging a program and find that the value of a variable is `None` when it shouldn't be, a function without a return statement is a likely cause.\n",
    "\n",
    "**Dysfunction #4:** Ignoring the return value.  Finally, here's a version where the function is correct, but the way it's used is not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def carrying_capacity(system):\n",
    "    K = -system.alpha / system.beta\n",
    "    return K\n",
    "    \n",
    "sys2 = System(alpha=0.025, beta=-0.0018)\n",
    "carrying_capacity(sys2)\n",
    "\n",
    "# print(K)     This line won't work because K only exists inside the function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, `carrying_capacity` runs and returns `K`, but the return value is dropped.\n",
    "\n",
    "When you call a function that returns a value, you should do something with the result.  Often you assign it to a variable, as in the previous examples, but you can also use it as part of an expression.\n",
    "\n",
    "For example, you could eliminate the temporary variable `pop` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(carrying_capacity(sys1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or if you had more than one system, you could compute the total carrying capacity like this:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "total = carrying_capacity(sys1) + carrying_capacity(sys2)\n",
    "total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:** In the book, I present a different way to parameterize the quadratic model:\n",
    "\n",
    "$ \\Delta p = r p (1 - p / K) $\n",
    "\n",
    "where $r=\\alpha$ and $K=-\\alpha/\\beta$.  Write a version of `update_func` that implements this version of the model.  Test it by computing the values of `r` and `K` that correspond to `alpha=0.025, beta=-0.0018`, and confirm that you get the same results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
